Heat conduction in a confined solid strip: Response to external strain 
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We study heat conduction in a system of hard disks confined to a narrow two dimensional channel. 
The system is initially in a high density solid-like phase. We study, through nonequilibrium molecular 
dynamics simulations, the dependence of the heat current on an externally applied elongational 
strain. The strain leads to deformation and failure of the solid and we find that the changes in 
internal structure can lead to very sharp changes in the heat current. A simple free-volume type 
calculation of the heat current in a finite hard-disk system is proposed. This reproduces some 
qualitative features of the current-strain graph for small strains. 
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PACS numbers: 62.20.Mk, 64.70.Dv, 64.60.Ak, 82.70.Dd 



I. INTRODUCTION 

In a recent study pj it was observed that the properties 
of a solid that is confined in a narrow channel can change 
drastically for small changes in applied external strain. 
This was related to structural changes at the microscopic 
level such as a change in the number of layers of atoms in 
the confining direction. These effects occur basically as 
a result of the small ( few atomic layers in one direction 
) dimensions of the system considered and confinement 
along some direction. A similar layering transition, in 
which the number of smectic layers in a confined liquid 
changes in discrete steps with increase in the wall-to- 
wall separation, was noted in Q, Q . Both Q, look at 
equilibrium properties while [j| looks at changes in the 
dynamical properties. An interesting question is, how 
are transport properties, such as electrical and thermal 
conductivity, affected for these nanoscale systems under 
strain? This question is also important to address in 
view of the current interest in the properties of nanosys- 
tems, both from the point of view of fundamentals and 
applications 0,11, El- 

In this paper we consider the effect of strain on the heat 
current across a two-dimensional (2D) "solid" formed by 
a few layers of interacting atoms confined in a long nar- 
row channel. We note here that, in the thermodynamic 
limit it is expected that there can be no true solid phase 
in this quasi-one-dimensional system. However for a long 
but finite channel, which is our interest here, and at a 
high packing fraction the fluctuations are small and the 
system behaves like a solid. We will use the word "solid" 
in this sense. 

In Ref. Q the anomalous failure, under strain, of a 
narrow strip of a 2D solid formed by hard disks confined 
within hard walls [ see Fig.^] wa s studied. Sharp jumps 
in the stress vs strain plots were observed. These were 



related to structural changes in the system which under- 
went transitions from solid-to-smectic-to-modulated liq- 
uid phases [jj, [jj . In the present paper we study changes 
in the thermal conductance of this system as it undergoes 
elastic deformation and failure through a layering tran- 
sition caused by external elongational strains applied in 
different directions. 

The calculation of heat conductivity in a many body 
system is a difficult problem. The Kubo formula and 
Boltzmann kinetic theory provide formal expressions for 
the thermal conductivity. In practice these are usually 
difficult to evaluate without making drastic approxima- 
tions. More importantly a large number of recent stud- 
ies |$, |2l [lfj, [llj indicate that the heat conductivity of 
low-dimensional systems infact diverge. It is then more 
sensible to calculate directly the heat current or the con- 
ductance of the system rather than the heat conductivity. 
In this paper we propose a simple-minded calculation of 
the heat current which can be expected to be good for a 
hard disk (or hard spheres in the three dimensional case) 
system in the solid phase. This reproduces some qualita- 
tive features of the simulations and gives values for the 
current which are of the correct order of magnitude. 

The organization of the paper is as follows. In Sec. |I1J 
we explain the model and present the results from simu- 
lations. In Sec. (|III|I we derive a simple formula for heat 
current in a hard-sphere system and evaluate it approxi- 
mately. We conclude with some discussions in Sec. (IIV(I . 



II. RESULTS FROM SIMULATIONS 

We consider a 2D system of hard disks of diameter d 
and mass m which interact with each other through elas- 
tic collisions. The particles are confined within a narrow 
hard structureless channel [see Fig.^. The hard walls of 
the channel are located at y = and y = L y and we take 
periodic boundary conditions in the a;— direction. The 
length of the channel along the x— direction is L x and the 
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area is A = L x x Ly. The confining walls are maintained 
at two different temperatures ( T2 at y = and T± at 
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FIG. 1: A solid with a triangular lattice structure formed by 
hard disks confined between two structureless walls at y = 
and y = L y . The walls are maintained at two different tem- 
peratures. The lattice parameters of the unstrained solid are 
denoted by a° x and a y . Elongational strains can be imposed 
by rescaling distances either in the x or y directions and the 
lattice parameters change to a x and a y . 



y = Ly ) so that the temperature difference AT = T2—T1 
gives rise to a heat current in the y-direction. Initially 
we start with channel dimensions L x and L® such that 
the system is in a phase corresponding to a unstrained 
solid with a triangular lattice structure. We then study 
the heat current in this system when it is strained (a) 
along the x— direction and (b) along the y— direction. 

We perform an event-driven collision time dynamics 
[l2| simulation of the hard disk system. The upper 
and lower walls are maintained at temperatures T% = 1 
and T2 = 2 (in arbitrary units) respectively by impos- 
ing Maxwell boundary condition Q at the two confining 
walls. This means that whenever a hard disk collides 
with either the lower or the upper wall it gets reflected 
back into the system with a velocity chosen from the dis- 
tribution 

where TV is the temperature (Ti or T2) of the wall on 
which the collision occurs. During each collision energy is 
exchanged between the system and the bath. Thus in our 
molecular dynamics simulation, the average heat current 
flowing through the system can be found easily by com- 
puting the net heat loss from the system to the two baths 
(say Qi and Q2 respectively) during a large time interval 
t. The steady state heat current from lower to upper 
bath is given by (I) = lim^oo Qi/r = - linv^oo Q2/T. 
In the steady state the heat current (the heat flux den- 
sity integrated over x) is independent of y. This is a 
requirement coming from current conservation. However 
if the system has inhomogeneities then the flux density 
itself can have a spatial dependence and in general we 
can have j = j(x,y). In our simulations we have also 
looked at j(x, 0) and j(x, L y ). 



FIG. 2: (Color online) Plots obtained by superposition of 500 
steady state configurations of a portion of 40 x 10 system taken 
at equal time intervals. Starting from r\ — 0.85 imposition of 
strains (a) e xx = 0.1, (b) e xx — 0.15 gives rise to these struc- 
tures. The colors code local density of points from red/dark 
(high) to blue/light (low). In (a) one can see a 9-layered struc- 
ture nucleated within a 10-layered solid. The corresponding 
structure factor identifies this to be a smectic P. In (b) the 
whole system has transformed into a 9-layered smectic. 



Note that the relevant scales in the problem are: ksT 
for energy, d for length and r s = \J md 2 /ksT for time. 
We start from a solid commensurate with its wall to wall 
separation and follow two different straining protocols. 
In case (a) we strain the solid by rescaling the length in 
the ^-direction and the imposed external strain is e xx = 
(L x — L x )/L x . In case (b) we rescale the length along the 
y-direction and the imposed strain is e yy — (L y — Ly)/Ly. 

The only thermodynamically relevant variable for a 
hard disk system is the packing fraction rj = ttNcP/AA. 
For a close packed solid with periodic boundary con- 
dition this value is about i] c = 0.9069. On the other 
hand for a confined solid having N y number of layers 
i] c = TrN y /(2V3{Ny - 1) + 4) and for a 10- layered solid 
7] c = 0.893. In our simulations we consider initial val- 
ues of rj for the solid to be close to t] c . The channel 
is "mesoscopic" in the sense that it has a small width 
with N y — 10 layers of disks in the y— direction (in the 
initially unstrained solid). In the x— direction the sys- 
tem can be big and we consider N x = 20, 40, 80, 160 
number of disks in the £— direction. In collision time dy- 
namics we perform 10 5 collisions per particle to reach the 
steady state and collect data over another 10 5 collisions 
per particle. All the currents calculated in this study are 
accurate within error bars which are less than 3% of the 
average current. 

Let us briefly mention some of the equilibrium results 
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FIG. 3: (Color online) Plot of j y vs e xx for different lengths of 
the channel. In this graph as well as in all other view graphs 
we plot j y in units of fcsT/r s d. The width of the channel is 
N y — 10 layers. Starting packing fraction is rj = 0.85. The 
solid line shows the theoretical prediction of dependence of 
the heat current on strain [see Sec. HIIIl l, 



for the stress-strain behavior obtained in Rcf. Q. As 
the strain e xx is imposed, the perfectly triangular solid 
shows rectangular distortion along with a linear response 
in strain vs stress behavior. Above a critical strain (e xx fa 
0.1) one finds that smectic bands having a lesser number 
of layer nucleate within the solid [this can also be seen 
in Fig. (2a) obtained from a nonequilibrium simulation] . 
This smectic is liquid-like in the x-direction (parallel to 
the walls) and has solid-like density modulation order 
in the y-direction (perpendicular to the walls). With 
further increase in strain, the size of the smectic region 
increases and ultimately the whole system goes over to 
the smectic phase at e xx fa 0.15 [Fig. (2b)]. At even 
higher strains the smectic melts to a modulated liquid 
0, 0] • The corresponding structure factor shows typical 
liquid like ring pattern superimposed with smectic like 
density modulation peaks. This layering transition is an 
effect of finite size in the confining direction. Similar 
phase behavior has been observed in experiments on steel 
balls confined in quasi ID ^^]. We note that, to fit a N, 
layered triangular solid within a channel of width L y 
require 
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This enables us to define a fictitious number of layers 

L v — d 
X = 2^=— + 1 



V3a 

of triangular solid that can span the channel where a is 
the lattice parameter at any given density. The actual 
number of layers that are present in the strained solid 
is N y — I(x) where the function I(x) gives the integer 



part of x- F° r confined solids the free energy has minima 
at integer values of x an d maxima at half-integral values 
P, Q • The difference in free-energy between successive 
maxima and minima gradually decreases with increas- 
ing L y . Thereby the layering transition washes out for 
ni > 25 layered unstrained solidQ. Up to this number 
of layers, a triangular solid strip confined between two 
planar walls fails at a critical deviatoric strain e* d ~ 1/iV^ 
(where td = £xx — e yy)- Smaller strips fail at a larger 
deviatoric strain. 

We now present the heat conduction simulation results 
for the two cases of straining in x and y directions. 

(a) Strain in ^-direction - In Fig. [21 we plot the 
heat current density j y calculated at different values of 
the strain e xx . Starting from the triangular lattice con- 
figuration, we find that the heat current decreases lin- 
early with increase in strain. At about the critical strain 
e xx fa 0.1 we find that the heat current begins to fall at a 
faster rate. This is easy to understand physically. At the 
onset of critical strain, smectic bands, which have lesser 
number of particle layers, start nucleating (Fig.|2J). These 
regions are much less effective in transmitting heat than 
the solid phase and the heat current falls rapidly as the 
size of the smectic bands grow. At about the strain value 
e xx fa 0.15 the whole system is spanned by the smectic. 
Beyond this strain there is no appreciable change in the 
heat current. The solid line in Fig.[3|is an estimate from 
a simple analysis explained in Sec. i|III|) . 

In Fig. 0] we plot the local steady state heat cur- 
rent j y (x) for a system of 40 x 10 particles at a strain 
txx = 0.118 i.e. at a strain corresponding to the solid- 
smectic phase coexistence. At this same strain the num- 
ber of layers averaged over 10 3 configurations have been 
plotted. It clearly shows that the local heat current is 
smaller in regions with smaller number of layers. This is 
the reason behind getting a sharp drop in average heat 
current after the onset of phase coexistence. 

(b) Strain in y-direction - Next we consider the case 
where, again starting from the density rj = 0.85, we im- 
pose a strain along the y— direction. As shown in Fig. 
the heat current j y now has a completely different nature. 
The initial fall is much steeper and has a form different 
from the linear drop in Fig. [3] The approximate analytic 
curve is explained in Sec. (fllll) . At about e yy fa 0.1 we 
see a sharp and presumably discontinuous jump in the 
current. At this point the system goes over to a buck- 
led phase (Fig. [BJd) in which different parts of the solid 
(along x-direction) are displaced along the y-direction by 
small amounts so that the extra space between the walls 
is covered 0, 0, 0]. A further small strain induces 
a layering transition and the system breaks into two re- 
gions one of which is an N y = 11 layered solid and the 
other is a N y — 10 layered highly fluctuating smectic- 
like region. At even higher strains (e yy ~ 0.2) the whole 
system eventually melts to an N y = 11 layered smectic 
phase. The phase behavior of this system is interesting 
and will be discussed in detail elsewhere^^. Unlike the 
case where the applied strain is in the x-direction, in 
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FIG. 4: (Color online) x( x ) is the local number of layers 
averaged over 10 3 steady state configurations for a system of 
40 x 10 hard disks. A starting triangular lattice of r] = 0.85 
is strained to e xx = 0.118 and the data collected after steady 
state reached. Also shown is the local heat current j y (x). The 
regions having lower number of layers conduct less effectively. 



the present case the buckling-layering transition is very 
sharp. Even though the overall density has decreased, 
due to buckling and increase in number of layers in the 
conducting direction, there is an increase in the energy 
transferring collisions and hence the heat current. The 
plots in Fig. show the structural changes that occur in 
the system as one goes through the transition. 

We find in general that the heat current along any 
direction within the solid shows the same qualitative fea- 
tures as the stress component along the same direction. 
This can be seen in Fig where we have plotted j y vs 
e xx for two starting densities of solids rj = 0.85, 0.89. In 
the inset we show the corresponding — a yy vs e xx curves 
and see that they follow the same qualitative behavior as 
the heat current curves. The reason for this is that mi- 
croscopically they both originate from interparticle col- 
lisions. Infact the microscopic expressions for the total 
heat current [ see Eq. in Sec. (jlllfl ] is very similar 
to that for the stress tensor component, with an extra 
velocity factor. The stress tensor is given by: 
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where {xf,uf} refer to the a-th component of position 
and velocity of the i th particle, r?- = J2 a ( x ?j) 2 anc ^ ^ij) 
is the interparticle potential. For a hard disk system, 
d ^ r '^ can be replaced by — fcsT(5(ry — d). Also in equi- 
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FIG. 5: (Color online) Plot of j y vs e yy for different channel 
lengths. The channel width is N y = 10 layers. The starting 
packing fraction is r\ = 0.85. The jump in current occurs at 
the strain value where the number of layers in the y— direction 
increases by one and the system goes to a smectic phase. The 
solid line shows the theoretical prediction of dependence of 
the heat current on strain [ see Sec. 1)1110 ]. 



Using collision time simulation it is easier to evaluate the 
stress tensor in the following way. We can rewrite Eq. © 
as 

We use the fact that (...) can be replaced by a time 
average so that from Eq. we have 
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Now note that during a collision we have J dtf^ — Apy 

where Ajjy is the change in momentum of i th particle 
due to collision with j th particle. It can be shown that 



Apy = -(uijJy)ry where fy = r^/ry and u 



; (3) Uj and r,, U; are evaluated just before a collision. This 



change in momentum occurs for a single pair of particle 
during one collision event. To get the stress tensor we 
sum over all the collision events in the time interval r 
between all pairs of particles. Therefore for collision time 
dynamics we get the following expression for the stress 
tensor, 



Acr a p — —Nk B TS, 
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where ^2 T denotes a summation over all collisions in time 

T. 
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FIG. 6: (Color online) Plots obtained by superposition of 500 
steady state configurations of a portion of 40 x 10 system taken 
at equal time intervals. Starting from r\ — 0.85 imposition of 
strains (a) e yy = 0.05, (b) e yy = 0.1, (c) e yy = 0.12 gives rise 
to these structures. The colors code local density of points 
from red/dark (high) to blue/light (low), (a) Solid phase, (b) 
A mixture of 10-layered solid and a buckling phase, (c) An 
11-layered solid in contact with 10-layered smectic like region. 



III. ANALYSIS OF QUALITATIVE FEATURES 

We briefly outline a derivation of the expression for the 
heat flux. For the special case of a hard disk system this 
simplifies somewhat. We will show that starting from this 
expression and making rather simple minded approxima- 
tions we can explain some of the observed results for heat 
flux as a function of imposed external strain. 

We consider a system with a general Hamiltonian given 
by: 
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V(vi)] 
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where V(rj) is an onsite potential which also includes 
the wall. To define the heat current density we need to 
write a continuity equation of the form: de(r,t)/dt + 
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FIG. 7: (Color online) Plot of j y vs e xx for two different start- 
ing values of the packing fraction. O corresponds to a starting 
value of T] = 0.89 while + is for r\ = 0.85. In both the cases 
the initial solid size was 80 x 10. The inset shows correspond- 
ing plots of — G yv (in units of ksT/d 2 ) vs e xx . Notice that 
stress-strain curve has the same qualitative profile as the j y 
vs e TT curve. 



dj a ( r ,t)/dx a — 0. The local energy density is given by: 
e(r, t) = '^^S(r — r i )h i where 

i 

Taking a derivative with respect to time gives 

''' .J-J2^-ri)hiuf + J2^-ri)hi (6) 
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where j K = J2i S{r — TijhiUi is the convective part of the 
energy current. We will now try to write the remaining 
part given by W u as a divergence term. We have 



^2S(r - r,-)fti 
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where /j" — —d(j){rij) / dxf . Using the equation of motion 
mil? = -dV/dxf + J2^ t /g we get 
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With the identification W u — —dj^/dx a and using 
we finally get: 



\ e f w - <o n ^ - + «?) ( 9 ) 



where #(2;) is the Heaviside step function. This formula 
has a simple physical interpretation. First note that we 
need to sum over only those i for which xf > x a . Then 
the formula basically gives us the net rate at which work 
is done by particles on the left of x a on the particles on 
the right which is thus the rate at which energy flows 
from left to right. The other part, j£ , gives the energy 
flow as a result of physical motion of particles across x a . 
Let us look at the total current in the system. Integrating 
the current density over all space we get: 
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Including the convective part and taking an average over 
the steady state we finally get: 



(i a ) = (i£) + (%)=^(hiu?) 



4 Ai. 
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(11) 



We note that for a general phase space variable 
A({xi, Ui}) the average (A) is the time average 
linv^l/r) J T dtA({xi(t),Ui(t)}). 

Finding the energy current for a hard disk sys- 
tem: The energy current expression involves the veloci- 
ties of the colliding particles which change during a col- 
lision so we have to be careful. We use the following 
expression for (ijf): 
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Now if we integrate across a collision we see that 
j di(fjj.Uj) gives the change in kinetic energy of the i th 
particle during the collision while J dt(fji.\ij) gives the 
change in kinetic energy of the j th particle. Hence we 
get 
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where we have used the fact that for elastic collisions 
AKi = —AKj and J^t denotes a summation over all 
collisions, in the time interval r, between pairs {ij}. 
The time interval between successive collisions between 
j th and j particles is denoted by and the aver- 
age ( ... ) c in the last line denotes a collisional aver- 
age. Thus (Tij) c = lim T ^ 00 r/ATy (t), where iVjj(r) is 
the number of collisions between i th and j th particles 
in time r. For hard spheres the convective part of the 
current involves only the kinetic energy and is given by 
{ la ) = J2i ( ( mu i/2) u i )• Using these expressions we 
now try to obtain estimates of the heat current and its de- 
pendence on strain (near the close packed limit where the 
system looks like a solid with the structure of a strained 
triangular lattice). 

Near the close packed limit the convection current 
can be neglected and we focus only on the conductive 
part given by (I u ) — (J^ 7 ) (for conduction along the 
y— direction). At this point we assume local thermal equi- 
librium (LTE) which we prove from our simulation data 
at the end of this section. Assuming LTE we write the 
following approximate form for the energy change AKi 
during a collision: 

AKi = k B (T{ Vj ) - T(yi)) = - k B^yi 3 = Vij kB ^ T > 

where we have denoted x\ a ~ 2 ^ = and AT = Ti — T\. 
The temperature gradient has been assumed to be small 
and constant. Further we assume that in the close packed 
limit that we are considering, only nearest neighbor pairs 
{< ij >} contribute to the current in Eq. (|13fl and that 
they contribute equally. We then get the following ap- 
proximate form for the total current: 
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where r c is the average time between successive collisions 
between two particles while y\ is the mean square sepa- 
ration along the y— axis of the colliding particles. Finally, 
denoting the density of particles by p = N/A we get for 
the current density: 
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(I v ) w 3pk B AT yj 

•A Ly T c 



(15) 



For strains e xx and e yy in the x and y directions we have 
P = Po/[(l + e xx) (1 + e yy)]- We estimate yl and t c from 
a simple equilibrium free-volume theory, known as fixed 
neighbour free volume theory (FNFVT). In this picture 
we think of a single disk moving in a fixed cage formed 
by taking the average positions of its six nearest neigh- 
bor disks [see Fig. |H]. For different values of the strains 
we then evaluate the average values [y1\f v and [r c ]/ for 
the moving particle from FNFVT. We assume that the 
position of the center of the moving disk Pq(x, y), at the 
time of collision with any one of the six fixed disks, is uni- 
formly distributed on the boundary B of the free- volume. 
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FIG. 8: In our free- volume theory we assume that the outer 
six disks are fixed and the central disk moves within this cage 
of fixed particles. The curve in bold line shows the boundary 
B of the free volume. A point on this boundary is denoted by 
Po(x,y) while the centers of the six fixed disks are denoted 
by Pi(xi,yi) with i = 1,2. ..6. 



Hence [y 2 ]f v is easily calculated using the expression: 



[vl 



fv 



(16) 



where Bt is the part of the boundary B of the free volume 
when the middle disk is in contact with the i th fixed 
disk, ds is the infinitesimal length element on B while 
Lb is the total length of B. Let the unstrained lattice 
parameters be a? x , a y = \/3a®/2. Under strain we have 
a x = a x (l + e xx ) and a v — a°(l + e aa ). Using elementary 
geometry we can then evaluate [yl]f v from Eq. (|16() in 
terms of e xx , e yy and the unstrained lattice parameter 
a x . An exact calculation of [r c ]/„ is nontrivial. However 

we expect [r c ]f v — c V^/T 1 / 2 where Vf v is the "free 
volume" [see Fig. |H] and c is a constant factor of O(l) 
which we will use as a fitting parameter. The calculated 
values for [y 2 ]f v and [t c ]/„ (see the appendix) are shown 
in Fig. Also shown are their values obtained from an 
equilibrium simulation of a single disk moving inside the 
free volume cage. Thus we obtain the following estimate 
for the heat current: 
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We plot in Fig. and Fig. (JSJ the above estimate of 
[j y ]fv along with the results from simulations. We find 
that the overall features of the simulation are reproduced 
with c = 0.42. 

For small strains we find (see the appendix) 
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(71 +72£ — 73£ 2 ) where e stands 
Pi, P2, 7i ; 72,73 are all pos- 
itive constants that depend only on a®. For r\ = 0.85 we 
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FIG. 9: 



(Color online) Plots showing comparison of the ana- 



lytically calculated values of [j/ c ]/u (in units of d ) and [t c ]/„ 
(in units of r s ) with those obtained from a free volume sim- 
ulation of a single disk moving within the free volume cage. 
The free volume corresponds to a starting unstrained trian- 
gular lattice at r\ — 0.85 which is then strained along x or y 
directions. 



have a = 7.62, ft = 121.77, ft = 124.37 and 71 = 0.02, 
72 = 0.33, 73 = 1.125. From these small strain scaling 
forms we find that j y (£ yy ) always decreases with positive 
e yy and increases with negative or compressive e yy (note 
that we always consider starting configurations of a tri- 
angular solid of any density) . On the other hand the sign 
of the change in j y (e xx ) will depend on the relative mag- 
nitudes of a, ft and 7 . For starting density 77 = 0.85, 
jyi^xx) decreases both for positive and negative e xx . In 
Fig-EH we show the effect of compressive strains e xx and 
e yy on the heat current j y and compare the simulation 
results with the free volume theory. 

It is possible to calculate y 2 and r c directly from our 
nonequilibrium collision time dynamics simulation. The 
mean collision time r c is obtained by dividing the total 
simulation time by the total number of collisions per col- 
liding pair. Similarly y 2 is evaluated at every collision 
and we then obtain its average. Inserting these values of 
t c and (y 2 ) into the right hand side of Eq. (|Trj|) we get 
an estimate of the current as given by our theory (with- 
out making use of free- volume theory). In Fig. 1111 we 
compare this value of the current j y , for strain e = e xx , 
and compare it with the simulation results. The excel- 
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FIG. 10: (Color online) Plot showing effect on j y of negative 
strains applied in the x and y directions. The system is pre- 
pared initially in a triangular lattice at r\ — 0.85. Note that 
negative e xx reduces j y whereas negative e yy increases j y . 



lent agreement between the two indicates that our simple 
theory is quite accurate. 

We have also tested the assumptions of a linear tem- 
perature profile and the assumption of local thermal equi- 
librium (LTE) that we have used in our theory. In our 
simulations the local temperature is defined from the lo- 
cal kinetic energy density, i.e. ksT = (mu 2 /2). Local 
thermal equilibrium requires a close to Gaussian distri- 
bution of the local velocity with a width given by the 
same temperature. The assumption of LTE can thus be 
tested by looking at higher moments of the velocity, eval- 
uated locally. Thus we should have (u 4 ) = 8(fc_BT/m) 2 . 
From our simulation we find out (u 4 (y)) and kBT(y) as 
functions of the distance y from the cold to hot reservoir. 
The plot in Fig. ^| shows that the temperature profile 
is approximately linear and LTE is approximately valid. 
We use our theory only in the solid phase and in this case 
there is not much variation in the direction transverse to 
the direction of heat flow (a;-direction). 

Finally we have verified that heat conduction in the 
small confined lattice under small strains shows a linear 
response behaviour. This can be seen in Fig. 1131 where we 
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FIG. 11: (Color online) Comparison of simulation results for 
j y with the approximate formula in Eq. I|15[l where r c and 
y\ are also calculated directly from the same simulation. The 
results are for a 40 x 10 system with starting value of r\ — 0.85 
and strained along £— direction. 




FIG. 12: (Color online) Plot of temperature profile and fourth 
moment of velocity for a strained 40 x 10 lattice. The un- 
strained packing fraction was r\ — 0.85 and the system was 
strained to e xx = 0.0625. 



plot j y vs AT = T 2 — Ti for a 40 x 10 triangular lattice at 
r\ = 0.85. Note that, as mentioned in the introduction, 
the bulk thermal conductivity of a two-dimensional sys- 
tem is expected to be divergent and the linear response 
behaviour observed here is only relevant for a finite sys- 
tem and in a certain regime (solid under small strain). 
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FIG. 13: (Color online) Plot of j y vs AT = T 2 - Tj for a 
40 x 10 triangular lattice at n — 0.85. We see that the current 
increases linearly with the applied temperature difference. 

IV. DISCUSSION 

In this paper we have studied heat conduction in a 
two-dimensional solid formed from hard disks confined 
in a narrow structureless channel. The channel has a 
small width (~ 10 particle layers) and is long (~ 100 
particles). Thus our system is in the nanoscale regime. 
We have shown that structural changes that occur when 
this solid is strained can lead to sudden jumps in the 
heat current. From the system sizes that we have stud- 
ied it is not possible to conclude that these jumps will 
persist in the limit that the channel length becomes in- 
finite. However the finite size results are interesting and 
relevant since real nano-sized solids are small. We have 
also proposed a free volume theory type calculation of the 
heat current. While being heuristic it gives correct order 
of magnitude estimates and also reproduces qualitative 
trends in the current-strain graph. This simple approach 
should be useful in calculating the heat conductivity of a 
hard sphere solid in the high density limit. 



The property of large change of heat current could be 
utilized to make a system perform as a mechanically con- 
trolled switch of heat current. Similar results are also 
expected for the electrical conductance and this is shown 
to be true at least following one protocol of straining in 
Ref . ^3 . From this point of view it seems worthwhile to 
perform similar studies on transport in confined nano- 
systems in three dimensions and also with different in- 
terparticle interactions. 
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APPENDIX A: CALCULATION OF [yl] }v AND V fv 

Using free volume theory, as explained in the text, we 
get the following expressions for [y1]f v and Vf v : 

[vl]fv = \(0 + i/> + 4>) + \ (sin2fl + cos2^ - - - j 

Vf V = -2(0 + ip + (f>) - (sin20 + sin2V> -sin 2<j>) 

+ 4a y (a x — cos (Al) 

where, 9 = skf~ (a x /2), i\) = sin (a x /2 — cos</>) and 
4> = ta,TT 1 (2a y /a x ) - cos -1 (Z/2) with I = ^J(a x /2) 2 + a\ 

(all lengths are measured in units of d). For strain along 
x-direction we have a x = 0^.(1 + t xx ) and a y = a y (= 
"\/3a2/2) while for strain along y-direction we have a y = 
a y (1 + c yy ) and a x = a x . From the above expressions 
we can obtain Taylor expansions of [y^]f v and [t c ]/„ = 
cVji 2 /T l l 2 , about the zero strain value. These give the 
expressions for a, fti, 02, 71, 72, 73 used in the text. 
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